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Abstract. - Various experiments have found a boundary slip in hydrophobic microchannel 
flows, but a consistent understanding of the results is still lacking. While Molecular Dynamics 
(MD) simulations cannot reach the low shear rates and large system sizes of the experiments, 
it is often impossible to resolve the needed details with macroscopic approaches. We model 
the interaction between hydrophobic channel walls and a fluid by means of a multi-phase 
lattice Boltzmann model. Our mesoscopic approach overcomes the limitations of MD simu- 
lations and can reach the small flow velocities of known experiments. We reproduce results 
from experiments at small Knudsen numbers and other simulations, namely an increase of slip 
with increasing liquid-solid interactions, the slip being independent of the flow velocity, and a 
decreasing slip with increasing bulk pressure. Within our model we develop a semi-analytic 
approximation of the dependence of the slip on the pressure. 



During the last century it was widely assumed that the velocity of a Newtonian liquid at 
a surface is always identical to the velocity of the surface. However, in recent years well con- 
trolled experiments have shown a violation of the no-slip boundary condition in sub-micron 
sized geometries. Since then, experimental [1] and theoretical works [2], as well as computer 
simulations [3-9] have tried to improve our understanding of boundary slip. The complex 
behavior of a fluid close to a solid interface involves the interplay of many physical and chem- 
ical properties. These include the wettability of the solid, shear rate, pressure, surface charge, 
surface roughness, as well as impurities and dissolved gas. Since all those quantities have to 
be determined very precisely, it is not surprising that our understanding of the phenomenon is 
still unsatisfactory. Due to the large number of different parameters, a significant dispersion 
of the results can be observed for ostensibly similar systems [1], e.g. observed slip lengths vary 
between nanometres [10] and micrometers [11] and while some authors find a dependence of 
the slip on the flow velocity [12], others do not [11,13]. Most computer simulations apply 
Molecular Dynamics (MD) and report increasing slip with decreasing liquid density [6, 7] or 
liquid-solid interactions [8,14], while slip decreases with increasing pressure [4]. These simula- 
tions are usually limited to some tens of thousands of particles, lengths scales of nanometres 
and timescales of nanoseconds. Also, shear rates are orders of magnitude higher than in any 
experiment [1]. We overcome these limitations using the lattice Boltzmann (LB) algorithm - 
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a powerful method for simulating fluid dynamics [15]. Rather than tracking individual atoms 
and molecules, the dynamics of the single-particle distribution function rj of mesoscopic fluid 
packets is described. In contrast to MD simulations, this method is less computationally de- 
manding and allows to simulate experimentally accessible length and time scales. Our ansatz 
differs from other LB approaches where slip is introduced by generalizing no-slip bounce back 
boundary conditions to allow specular reflections with a given probability [3] or where the 
viscosity is modified due to local density variations [16]. In both cases, parameters deter- 
mining the properties at the boundaries are not easily mappable to experimentally available 
values. Our approach is based on Shan and Chen's multi-phase LB model [17]. Here, inter- 
actions between different species are modelled by mesoscopic forces between the phases. This 
naturally opens the way to introduce similar interactions between each fluid species and the 
channel walls, where the strength of the interaction is determined by the fluid densities, free 
coupling constants, and a wall interaction parameter which is treated in a similar manner as 
a local fluid density. The model allows the simulation of multi-phase flows along hydrophobic 
boundaries and is introduced in the following. However, in order to study the influence of 
hydrophobicity on the boundary slip and to demonstrate the basic properties of the model, we 
focus on single phase flow in this paper. Results of multi-phase simulations will be presented 
in a future work. A multi-phase LB system can be represented by a set of equations [18] 

7 ? f(x + c i ,t + l)-7 7 f(x,t) = Of, i = 0, 1, . . . , 6 , (1) 

where rjf (x, t) is the single-particle distribution function, indicating the amount of species a 
with velocity c^, at site x on a D-dimensional lattice of coordination number b (D3Q19 in 
our implementation), at time-step t. For the collision operator Qf we choose the Bhatnagar- 
Gross-Krook (BGK) form 

=-^(v?(x,t)-V? e9 (u a (x,t),r ] <*(x,t))) , (2) 

where t q is the mean collision time for component a and determines the fluid viscosity. The 
system relaxes to an equilibrium distribution rj" eq which can be derived imposing restriction on 
the microscopic processes, such as explicit mass and momentum conservation for each species 
[19]. jy a (x, t) = J2i Vi ( x j i s the fluid density and u a (x, f) is the macroscopic velocity of the 
fluid, defined as r] a (x, t)u a (x, t) = J2i V? ( x j *) c »- Interactions between different fluid species 
are introduced as a mean field body force between nearest neighbors [17]: 

F Qa (x,t) S -^(x,i)^.g aa ^V a (x',i)(x'-x) , (3) 

where ^ a (x,t) = (1 — e~ r ' ( x >*)/'ro) is the so-called effective mass with 770 being a reference 
density that is set to 1 in our case [17]. g &a is a force coupling constant, whose magnitude 
controls the strength of the interaction between component a and a. The dynamical effect of 
the force is realized in the BGK collision operator in Eq. J5J) by adding to the velocity u in 
the equilibrium distribution an increment <5u" = T a F aa /rj a . For the interaction of the fluid 
components with the channel walls we apply mid-grid bounce back boundary conditions [15] 
and assign interaction properties to the wall which are similar to those of an additional fluid 
species. I.e., we specify constant values for the force coupling constant g &a = 9waii,a and the 
density r] a = r/ waU at wall boundary nodes of the lattice. This results in a purely local force 
as given in Eq. between the flow and the boundaries. Even though one could argue that a 
single parameter to tune the fluid- wall interaction would be sufficient, we keep our approach as 
close as possible to the original idea of Shan and Chen in order to benefit from the experience 
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obtained from other works using the original model. Furthermore, the additional parameter 
allows more flexibility when simulating not only a single fluid, but a multi-phase system. The 
fluid-wall interaction can be linked to a contact angle between fluid droplets and solid walls as 
it is often used to quantitatively describe hydrophobic interactions [20] . Recently, Benzi et al. 
have shown how to compute the contact angle within the Shan-Chen model [21]. The same 
authors also developed an approach to model apparent slip which is related to ours, but instead 
of using only local fluid-solid interactions, they add an exponential decay of the interaction 
with distance from the wall [22]. We simulate pressure driven flow between two infinite planes 
(Poiseuille flow) , where pressure driven boundary conditions are implemented in a similar way 
as in most experiments: a fixed pressure is set at the channel inlet and an open boundary at 
the outlet. The outlet is realized by interpolating the particle distribution function at the end 
of the channel as given by rff (x, t) = 2rjf (x — 1, t) — r]f (x — 2, t) leading to a linear pressure 
gradient. Already in 1823, Navier proposed a boundary condition where the fluid velocity 
at a surface is proportional to the shear rate at the surface, i.e. v z (xq) = (3dv z (x)/dx at 
x — xq [23]. Following his hypothesis the velocity in flow direction (v z ) at position x between 
the planes is given by 

where 2h is the distance between the planes, and p the viscosity. In contrast to a no-slip 
formulation, the last term in Eq. linearly depends on the slip length (3. Since (3 is typically 
of the order of nanometers or micrometers, it can be neglected in macroscopic experiments. 
In order to obtain (3 from our data, we measure the pressure gradient dP/dz at the center of 
the channel and the velocity profile between the two planes at a fixed position z. /3 is then 
obtained by a least square fit with Eq. 0] 

Our simulation parameters are as follows: the lattice size is kept constant with the channel 
length (z direction) being 256 sites, the distance between the plates 2h being 60 sites (x 
direction). We approximate infinite planes by using a 16 sites wide channel with periodic 
boundaries in y direction. In order to assure a fully equilibrated system we simulate for at 
least 40000 time steps before measuring and assured our results being independent of the 
discretization level by comparing to simulations of 28 and 124 sites wide channels. Each data 
point in the figures below corresponds to about six hours simulation time on eight IBM Power 
4 1.7GHz CPUs. All units in this paper are in lattice units with the lattice constant c and 
timestep At set to 1 if not stated otherwise. 

The dependence of the slip length (3 on the interaction parameter g wa ii, a is studied for 
r} wall =l.Q and 5.0. The bulk pressure P = pc 2 s , where p is the fluid density and c s — 1/V3 
the speed of sound, is kept at P=0.11, while the flow velocity is set to y=0.033. As shown in 
Fig.^, we vary g wa ii, a from 0.06 to 0.22 and find a steady increase of /3 for increasing g wa u, a - 
As expected, the curve for 7^^'"^ =5.0 is growing substantially faster than for r) wall = 1.0. The 
maximum available (3 are at about 5.2 for <7 lua ;;,a=0.26 and 77""" —1.0. At these strong fluid- 
wall interactions, the force as given in Eq. |3| becomes very large and results in a large area of 
low fluid density close to the wall. Increasing the interaction even further results in numerical 
instabilities due to too steep density gradients. In order to study the dependence of the slip on 
other parameters, the coupling constant g wa ii, a is kept constant at 0.08 from now on. FigQ]) 
depicts the dependence of (3 on r) wal1 for different bulk pressures P=0.033, 0.1, and 0.3 and 
fixed flow velocity V = 3.5 • 10 -3 in the system. While all three graphs grow constantly with 
increasing rj wal1 , the one for P=0.033 grows the fastest demonstrating that absolute values 
for (3 are higher for lower pressure. 

We have measured the magnitude of the boundary slip over a very wide range of flow 
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Fig. 1 - Slip versus g wa u, a for different wall interactions rj w and constant P=0.11, V=0.033 (a). f3 
is steadily increasing with increasing g wa ii,a and achievable slip lengths are higher for a larger r/ wal1 . 
Fig. b) shows (3 versus interaction parameter r] wal1 for various bulk pressures and fixed V = 3.5- 10 -3 . 
For lower pressure, larger values of f3 are measured. 



velocities V from 1 • 1CT 4 to 3 • 10" 2 for wall interactions Tf aU =0.0, 0.5, 1.0, and 2.0. V 
is measured at the center of the channel and given on a logarithmic scale in Fig. [21 For 
r jwaii_Q q we n0 £ £ n( ^ an y boundary slip confirming that our method properly reproduces 
no slip behavior in the interaction free case. With increasing wall interactions, we achieve 
an increase of the magnitude of /3 to up to ~1.1 for rf" =2.0. We are not able to find 
any velocity dependence of j3, but find constant slip for fixed fluid-wall interactions, which 
is consistent with many experiments [13,24]. The fluctuations of the data for very low flow 
velocities are due to numerical uncertainties of the fit at very low curvature of the parabolic 
velocity profile. For V > 0.01 we find a slight deviation of /3 from the constant measurements. 
This is due to a small variation of the bulk pressure from P=0.097 for V = 1- 10~ 4 to P=0.106 
for V = 0.03 that cannot easily be avoided for technical reasons. We have checked for a few 
data points that j3 stays constant if P can be kept at exactly fixed values, too. The slip length 
being independent of the flow velocity is consistent with many experiments and computer 
simulations, like the MD simulations of Cottin-Bizonne et al. [25] and the experiments of 
Cheng et al. [13] and Baudry et al. [26]. We speculate that an increase of (3 with increasing 
flow velocity as measured by some experiments [12] is due to surface roughness of the channel 
boundaries or other nonlinear effects. Since our model is not able to treat roughness on an 
atomic scale, we do not expect to conform with those results. MD simulations which find 
a non-constant value for j3 operate at very high shear rates which are orders of magnitude 
higher than what can be obtained by our approach [9] . 

Computing the exact slip in dependence of the interaction parameters from first principle 
analytically is a very hard or even impossible task since our interaction as given in Eq. 
modifies the equilibrium distribution in the BGK operator. Therefore, we present a semi- 
analytic approximation which utilizes the common two-layer model. Here, it is assumed that 
a thin fluid layer with thickness 6 and different viscosity as the bulk fluid exists near the 
channel walls. As calculated by various authors [2], within this model the slip length can be 
computed as (3 = {pbuik/p^i — !)<>> where fibuik is the viscosity of the bulk fluid, and p\ the 
viscosity close to the wall. Since the dynamic viscosity is given by the kinematic viscosity 
times the fluid density, jj, = pv = p(2r Q — l)/6 [15], we write j3 = (pbulk/pi — !)<$• Pbuik can be 
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Fig. 2 - Slip length (3 versus flow velocity V for different wall interactions r/ waU . While we do not 
find any slip for rf aU =0.0, (3 increases with increasing 77 . We vary the flow velocity from 1 • 10 -4 
to 0.03 and find constant values for /3 independent of V (within numerical accuracy). 

measured in the channel's center and p\ at the first lattice site next to the wall. Fig.Olshows 
the dependence of p\ on g wa ii, a for r^ — 1.0, 5.0, P=0.11, and V=0.033. Since pi cannot 
easily be computed analytically, we postulate an interaction term that depends on the bulk 
density and the fluid-wall interaction as well as a free fit parameter k, 



1 = kF 



wall, a 



(■K,t)/p bu lk{x,t) 



(5) 



and fit p\ with an exponential function p\ = pb«ifc(x, t) exp(— X s ). With only a single value for 
k we are able to utilize this equation to fit p\ for all our simulation parameters, k is found 
to be 8.35 for our data. The lines in Fig. |31 illustrate the good quality of our approximation. 
A similar approach is applied to model the thickness of the layer at the wall which strongly 
depends on the fluid- wall interaction and bulk density. Here, we set 5 = exp(X). As a result, 
[3 can be estimated by (3 — (exp(X) — 1) exp(X). The semi- analytic approximation is used to fit 
the dependence of the slip length (3 on the bulk pressure P. Fig. 0] shows the simulation data 




0.15 



0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 
Smill.a 



Fig. 3 - The fluid density close to the channel walls pi over g wa ii,c, is given for rf 
The lines correspond to a fit by our semi-analytic approximation. 



=1.0, 5.0 (symbols). 
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Fig. 4 - Slip length /3 versus bulk pressure P for jf" = 0.5,1.0,2.0 (symbols). (3 increases with 
increasing fluid- wall interactions, but decreases with increasing P. This dependence can be described 
by a semi-analytic equation (lines) which agrees well for small fluid- wall interactions and qualitatively 
reproduces our data for strong fluid-wall interactions. 



(symbols) and the approximation (lines) for wall interactions rf =0.5, 1.0, and 2.0. The 
bulk pressure is varied from 0.03 to 0.33. We find a decrease of (3 with increasing pressure 
P. An increase of rj waU leads to an increasing slope of the curves and to higher absolute 
values for (3. Furthermore, we find a decrease of the slip with increasing bulk pressure. These 
results qualitatively agree with MD simulations [4,8]. Even with a single value for the fit 
parameter k, the semi-analytic description of (3 agrees very well for low fluid-wall interactions. 
For strong interactions (rj wall =2.0), the fit qualitatively reproduces the behavior of the slip 
length. Higher order terms in the exponential ansatz for 5 are needed for a better agreement. 

To demonstrate that our approach is able to achieve experimentally available length and 
time scales, we scale our simulations to the experimental setup of Tretheway and Meinhart 
[11]. They use a 30/im high and 300/xm wide microchannel with typical flow velocities of 
V = 10 mm/s. For water, they measure a slip length of 0.92/im. The Reynolds number 
Re = 2hV/v in their experiment is ~0.3. To reproduce the observed slip, we set gw;,a=0.16 
and r] waU = 1.0 (see Fig. QJi). We arc able to cover a wide range of flow velocities, i.e. 
for the setup given above, velocities can range from as low as 1 x 10~ 4 to as high as 0.05 
corresponding to Re between 0.038 and 19. The Knudsen number is given by Kn = v/(c s 2h) 
which corresponds to 4.8 x 10~ 3 for the simulations presented here. At these low Kn, the 
hydrodynamic approach is well valid. However, it has been shown that the LB method can 
be applied for Kn much larger than 1 if one uses modified boundary conditions [27]. Our 
mesoscopic force is expected to be able to properly describe fluid-wall interactions in such 
systems as well. 

In conclusion, we presented a new approach to investigate boundary slip in hydrophobic 
microchannels by means of a multi-phase LB model. In contrast to MD simulations, our 
model is able to reach the length and time scales of typical experiments and is applicable for 
a wide range of realistic flow velocities. We qualitatively reproduced the dependence of slip 
on the hydrophobicity of the channel walls and found constant slip for varying flow velocities. 
The decrease of the slip with increasing pressure can be approximated by a semi-analytic 
approach. Our results are consistent with MD simulations [4,8,25] and experiments [11]. 
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